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I . INTRODUCTION 


Fragment  penetration  of  fuel  cells  is  a traditional  damage  and  kill 
mechanism  of  the  Army's  motorized  vehicles  and  fuel  storage  containers. 

In  order  to  provide  information  which  will  assist  in  understanding  the 
physical  phenomena  involved,  it  was  decided  that  a new  approach  for 
simulating  fluid  flow  may  be  of  value.  Consequently,  a completely  general 
formulation  was  conceived  to  serve  as  the  basis  for  the  development  of  a 
computer  code . ^ 

The  purpose  of  this  report  is  to  document  the  basic  ideas  given  in 
the  original  notes  and  to  present  an  expansion  of  them  in  a form  suitable 
for  use  as  a guide  in  the  construction  of  a code.  It  is  not  the  inten- 
tion of  this  report  to  discuss  all  of  the  physics  which  would  be 
involved  in  such  an  undertaking.  Presented  here  is  a mathematical 
description  of  the  procedure  for  simulating  fluid  flow  by  combining 
Lagrangian  cell  motion  and  the  computation  of  distributions  of  physical 
quantities  over  cells  based  on  a knowledge  of  moments  over  subregions 
of  the  cell. 


II.  THE  BASIS  OF  THE  PROPOSED  NUMERICAL  TECHNIQUE  FOR 
TRANSPORTING  PHYSICAL  QUANTITIES  IN  FLUID  FLOW 

The  numerical  technique  for  simulating  fluid  flow  transport  consists 
of  two  separate  parts;  the  sum  of  which  will  yield  distributions  of  the 
physical  parameters  across  fixed  cells  as  a function  of  time.  While  the 
discussion  is  presented  in  terms  of  rectilinear  coordinates,  the  same 
theory  is  applicable  in  cylindrical  or  spherical  coordinates  and, in  fact, 
the  code  should  be  constructed  to  accomodate  all  three.  To  facilitate 
understanding,  the  discussion  will  concentrate  on  one  part  of  the  numer- 
ical technique  at  a time. 

Figure  II-l  consists  of  a diagram  which  portrays  cells  assumed  to  be 
formed  by  a fixed  rectangular  coordinate  system  x,  y.  The  two  dimen- 
sional assumption  implied  by  the  figure  is  not  necessary  since  the 
arguments  are  equally  valid  in  one,  two,  or  three  dimensions.  In  each 
cell  an  arbitrary  physical  parameter  is  represented  by  the  symbol  Oj. 

The  subscript J is  to  denote  that  the  distribution  of  a corresponding  to 
each  J cell,  can  be  different.  This  situation  is  assumed  to  exist  at 
time  tg,  which  would  be  the  instant  corresponding  to  that  at  the  beginn- 
ing of  a cycle  of  computation. 

In  each  of  the  four  cells  depicted  in  Figure  II-l,  the  fluid  has 
associated  with  it  a local  velocity  distribution  of  the  form  presented 
in  Figure  II-2.  As  indicated,  the  two  components,  and  ^ , are  assumed 


^Rogers,  Joel,  unpublished  progpess  veports  submitted  to  the  Ballistic 
Research  Laboratory,  March  9,  1970  and  January  19,  1971. 


to  be  linear  and  in  general  are  represented  by  the  following  expres- 


= i (a^  y . a^) 

= 1 (a^  X ^ ^4^ 


(II-l) 


where  i and  j are  unit  directional  vectors  and  a , a-,  a , and  a are 
coefficients . 

Treating  the  fluid  in  the  cells  of  Figure  II-l  as  Lagrangian  cells 
which  move  as  bodies  of  fluid  according  to  the  fluid's  local  velocity, 
the  new  positions  of  these  Lagrangian  cells  at  tg  + At  can  be  obtained 
and  are  represented  by  the  dashed  lines  in  Figure  II-3.  The  symbol  At 
represents  a finite  interval  of  time.  The  shaded  regions,  denoted  by 
SRj,  constitute  contributions  of  the  moments  of  physical  parameters 
from  the  four  Lagrangian  cells  to  one  cell  in  the  fixed  coordinate 

system.  The  immediate  task  is  to  explain  how  to  obtain  the  distribu- 

tions of  the  physical  parameters  based  on  these  four  contributions. 

The  basis  for  the  procedure  for  computing  the  distribution  of  the 
physical  parameter  is  presented  in  Appendix  A.  There  is  shown  that  if 
certain  moments  of  the  physical  parameter  are  available  for  the  sub- 

regions,  then  a new  distribution  based  on  these  moments  can  be  calcu- 

lated for  that  parameter  across  the  total  region.  To  accomplish  this, 
it  is  necessary  to  resolve  the  question  as  to  which  moments  are  required 
and  how  are  the  moments  obtained?  For  the  present,  we  may  hold  the 
latter  portion  of  the  question  in  abeyance  and  assume  that  the  moments 
are  available.  The  decision,  relative  to  the  first  part  of  the  question 
above,  depends  on  the  form  of  the  distribution  of  the  physical  parameter 
desired  as  we  will  demonstrate  by  the  following  example. 

First  of  all,  we  may  arbitrarily  assume  that  the  distribution  of 
the  physical  parameter  a,  in  Figure  II-3,  is  to  be  of  a linear  form  as 
follows: 


a = bjX  + b2y  + b^ 


(11-2) 


where  a is  a function  of  the  spatial  parameters  x and  y and  bj^,  b2,  and 
b^  are  coefficients  to  be  evaluated.*  In  order  to  evaluate  the  three 
coefficients,  it  is  sufficient  to  generate  three  equations  containing 
these  coefficients  and  to  solve  them  simultaneously.  To  this  end,  we 
multiply  through  Equation  I I -2  by  a function  f which  is  required  to 


"The  linear  assumption  assumed  here  is  not  necessary  in  that  a higher 
order  polynomial  may  be  utilized.  In  any  case  the  same  procedure  is 
to  be  followed. 


take  on  values  of  1, 


and 


X,  and  y.  That  is: 

fa=bjfx+b2fy+b2 

a = bi  X -K  b2  y + bj 
X a = b^  x^  + b2  X y + bj  X 
2 

ya  = bj^xy  + b2y  +b^y. 


Cn-3) 


(11-4) 


Integrating  the  terms  of  Equations  I I -4  with  respect  to  x and  y 
yields  the  following  equations: 


y*a  dA  = bj  y*x  dA  + b2  dA  * bj  JdA 

'rc  TC  TC  TC 

J'x  a dA  = dk  + b^  Jx  y dk  * b^  Jx  dA 

rc  TC  TC  TC 

fy  a dA  = bj^  x y dA  + b2  J" y^  dA  + b^  ^y  dA 


TC 


where  dA  equals  dx  dy  and  the  symbol  TC  is  intended  to  indicate  that  the 
integrations  are  to  be  performed  over  limits  which  span  the  entire  cell 
(from  x„i„  and  y to  x^  and  y respectively  i"  Figure  II-l).  For 
other  forms  of  the  distribution,  the  function  f would  be  set  equal  to 
other  values,  but  in  general; 


f = X y 


(II-6) 


where  kj  and  k2  are  equal  to  0,  1,  2,  ....  For  combinations  of  these 
values  for  k]^  and  k2,  the  values  of  f are  1,  x,  y,  x y,  x^y,  xy^,  x y^ , 
and  etc. 


The  integrals  on  the  right  side  of  Equations  I I -5  depend  on  the 
limits  defined  by  the  coordinate  values  of  the  corner  points  of  the  cell 
and  are  easily  evaluated.  We  refer  to  them  as  Geometric  Moments  because 
their  values  depend  on  the  spatial  parameters  only,  and  for  convenience. 


*The  proposed  distributions  for  density ^ momentum^  and  kinetic  energy 
and  the  required  moments  are  presented  in  detail  in  Appendix  C. 
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wo  roiirosont  them  by  the  symbol  G as  follows: 

S ■ /" 

TC 

G2  = /y 

TC 

S = 

re 

(11-7) 

^4  = f^y  dA 

TC 

G5  = p dA 
TC 

G6  =/»'<«■ 

TC 


Substitution  of  Equations  II-7  into  Equations  II-5  yields  the 
following  equations: 


dA  = bi  G3  + b2  G2  + b^Gj 

a dA  = b^  G3  + b^G^  + b^G^ 

a dA  = b G + b G,  + b G . 
14  Z O 6 Z 


(11-8) 


The  integrals  on  the  left  are  equal  to  the  sums  of  the  integrals 
over  the  various  subregions.  Consequently,  the  following  equations  can 
be  written: 


14 
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■r-w 


/“  “ ■ / “l 


dA  + / dA  . / a3  dA  . / 

“4 


^ ^ ct  dA  = X dA  + J'  X 02  dA  + x dA  + ^ x 


(11-9) 


J" f a dk  - J"  y Oj  dA  + J'  y 0.^  dA  + J"  y dA  + J"  y dA 


where  SRj  indicates  that  the  limits  of  integrations  extend  over  the 
various  subregions. 

We  obtain  the  following  equations  after  combining  Equations  I I -8 
and  Equations  I 1-9: 


E/ 


ttj  dA  = bj  G3  + b2  G2  + b^  G^ 


T .'•n 


dA  = b,  Gp  + b_  G.  + b_  G_ 
1 5 2 4 3 3 


(11-10) 


E/^  aj  <iA  = bj  • ‘■j  <=2  • 


The  remaining  task  consist  of  evaluating  the  moments  of  the  physical 
parameter  over  the  various  subregions  as  indicated  by  the  integrals  on 
the  left  side  of  Equations  II-IO  and  the  shaded  regions  in  Figure  II-3. 

The  limits  of  these  subregions  can  be  determined  by  the  intersections 
of  the  boundaries  of  the  Lagrangian  cell  and  the  fixed  coordinate  system's 
cells.  However,  the  distributions  of  Uj  are  unknown  since  changes  in 
otj  have  occurred  over  the  time  step  At,  and  thus  the  integrals  cannot  be 
evaluated  directly. 

In  Figure  II-4A,  a Lagrangian  cell  of  fluid  is  shown  at  t^  + At, 
where  the  subregions  are  formed  by  the  intersections  of  the  boundaries. 
Since  the  motion  of  the  cell  is  calculated  using  the  local  fluid  velocity 
and  the  time  step  At,  the  limits  of  these  subregions  can  be  converted  by 
the  reverse  process  to  limits  of  corresponding  subregions  for  the  cell 
located  at  tg,  as  shown  in  Figure  II-4B.  The  distribution  of  a is  known 
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at  tg,  thus  if  the  integrals  for  calculating  moments  could  be  evaluated 
over  the  limits  at  tg.  such  that  the  results  are  equivalent  to  those  at 
tg  + At,  then  the  difficulty  would  be  resolved. 

In  order  to  permit  integrations  at  tg,  the  function  f(x,y,t)  is 
defined  at  tg  + At  as  follows: 

k k 

f(x,  y,  t^  + At)  = X ^ y ^ (II-ll) 

and  is  required  to  satisfy  the  following  differential  equation: 


. V f = 0 . (11-12) 

This  differential  equation  is  to  be  interpreted  as  saying  that  f is  a 
constant  along  particle  paths.  Thus,  to  first  order  accuracy,  we  have 


f(x,  y,  t^)  = f (X  + At,  y + At,  tg  + At)  (11-13) 


where  x,  y are  the  global  coordinates  of  the  fixed  coordinate  system. 


In  Appendix  B,  the  Equation  11-12  and  the  conservation  equations 
of  change  are  utilized  in  deriving  the  appropriate  equations  of  change 
needed  to  compute  moments  of  the  physical  parameters  over  the  subregion 


III.  A COMPUTATIONAL  VERIFICATION  OF  THE  METHOD 

A satisfactory  verification  of  the  method  proposed  in  this  report 
demands  a comparison  with  computations  of  complicated  problems  solved 
with  existing  codes  using  already  "proven"  methods.  This  is  currently 
impossible  since  the  program,  utilizing  this  proposed  method,  has  not 
been  constructed.  Consequently,  the  best  one  can  do  at  this  point  is  to 
solve  a simple  problem  for  the  purpose  of  demonstrating  that  the  pro- 
posed method  for  simulating  fluid  flow  is  feasible.  In  addition,  such 
a calculation  may  be  helpful  in  providing  a means  whereby  the  method 
can  be  better  understood.  In  any  case,  this  section  is  intended  to 
verify  the  method  to  a limited  extent  and  should  be  approached  from  that 
point  of  view. 

The  problem  chosen  is  the  determination  of  the  steady  velocity 
distribution  for  flow  between  two  semi-infinite  parallel  flat  plates, 
where  the  driving  mechanism  is  the  frictional  forces  due  to  the  constant 
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motion  of  one  of  the  plates.  Diagrams,  describing  the  temporal  aspect 
of  the  problem,  are  presented  in  Figure  III-l.  In  the  diagram  at  the 
top  of  Figure  III-l,  the  velocity  is  zero  throughout  the  fluid  for  times 
prior  to  some  arbitrary  instant  of  time  tQ.  At  tQ,  the  lower  plate  is 
set  in  motion  suddenly  at  some  velocity  (in  this  case  the  value  0.8  is 
assumed)  and  due  to  the  frictional  force  between  the  plate  and  the  fluid, 
the  fluid  is  dragged  along  with  the  plate.  It  is  assumed  that  no 
slippage  exist  between  the  surfaces  of  the  fluid  and  the  plates;  thus 
that  layer  of  fluid  adjacent  to  the  surface  of  the  plates  and  the  plates 
have  the  same  velocity.  For  all  time  after  tQ,  the  lower  plate  is 
forced  to  retain  its  initial  velocity  and  momentum  is  continuously 
transferred  in  the  y direction  from  the  moving  plate  to  the  fluid  and 
through  the  fluid  toward  the  upper  plate.  The  momentum  is  transferred 
from  one  layer  of  fluid  to  the  next  due  to  the  frictional  force  created 
by  the  effect  of  the  viscosity  of  the  fluid.  Since  the  upper  plate  is 
forced  to  remain  motionless,  the  velocity  of  the  layer  of  fluid  adjacent 
to  its  surface  remains  zero.  Eventually  the  steady  velocity  distribu- 
tion is  obtained,  as  shown  in  the  bottom  diagram  of  Figure  III-l.  The 
objective  of  the  calculation  is  to  determine  momentum  distributions  of 
the  fluid  from  tQ  until  the  steady  velocity  distribution  is  obtained. 

It  is  not  necessary  to  simulate  a real  fluid  to  make  our  point,  thus 
simple  values  are  used  for  the  parameters  involved  and  no  reference  to 
specific  units  are  made. 


The  fluid  motion  in  this  flow  problem  is  in  one  direction  only; 
that  being  in  a direction  parallel  to  the  plates.  The  form  of  the 
velocity  distribution  is  written  as  follows: 

v^  = t (a^  y + (III-l) 


where 

v^  = the  velocity  parallel  to  the  plates. 

y = the  coordinate  perpendicular  to  the  plates 
in  the  j direction. 

T = the  unit  vector  in  the  direction  parallel 
to  the  plates. 

^1’  ^2  ' coefficients. 

The  motion  of  the  Lagrangian  cell  is  demonstrated  in  Figures  III-2 
and  III-.3,  where  the  dashed  lines  represent  the  cells  formed  by  the 
fixed  coordinate  system  and  the  solid  lines  indicate  the  positions  of 
the  Lagrangian  cells  at  t + At.  Focussing  attention  on  the  T row  of 
cells  in  Figure  III-2,  we^see  that  the  cells  of  fluid  are  partitioned 
into  two  subregions  by  the  intersection  of  boundaries.  The  shaded 
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subregions  represent  the  contributions  of  fluid  from  Lagrangian  cells, 
which  were  initially  superimposed  on  (S-1,T)  and  (S,T)  fixed  cells  at 
tg,  to  the  fixed  cell  (S,T)  at  t + At.  Based  on  this  information,  the 
fixed  cells,  to  which  the  contributions  from  the  various  subregions  are 
to  be  summed,  are  identified.  The  motion  of  the  Lagrangian  cells  are 
computed  by  translating  the  comer  points  with  the  Equation  III-l  and  a 
time  step  At. 

_^In  this  problem,  at  any  point  in  time,  the  flow  does  not  vary  in 
the  i direction  parallel  to  the  flat  plates.  Consequently  the  velocity 
distribution  at  tg  + At  need  be  determined  in  only  one  column  of  cells 
(say  Column  S of  Figure  III-2).  However,  the  procedure  required  an 
additional  column  of  cells  (Column  S-1  in  Figure  III-2)  in  order  to 
obtain  the  contribution  of  flow  from  one  column  of  cells  into  the  next 
column.  In  this  case  the  contributions  are  represented  by  the  subregions 
identified  as  SRg.j .j  and  SR5  .p.  At  the  end  of  a cycle  of  computations 
both  columns  of  cells  are  given  the  same  velocity  distribution  as  that 
calculated  in  order  to  start  the  next  cycle  of  computations. 

For  this  problem,  only  momentum  need  be  considered  to  obtain  the 
required  results.  And  for  simplification,  we  assume  that  the  density 
distribution  throughout  the  flow  is  equal  to  1.  Therefore,  the  expres- 
sion for  the  momentum  is  identical  to  the  velocity  distribution  which 
is  expressed  as  Equation  III-l.  Thus,  the  two  moment  equations  needed 
to  obtain  the  velocity  or  momentum  distribution  across  the  fixed  cell 
at  tg  + At  are  written  as  follows: 

J dx  dy  + dx  dy  = aj  Jy  dx  dy  + /-  dy 


S-1,T 


y v^  dx  dy  + 


(111-2) 


dx  dy  = a^  f y dx  dy  + a^  / y dx  dy  . 


S-1,T 


To  evaluate  the  left  side  of  Equations  I I 1-2,  it  is  necessary  to 
revert  to  the  subregion's  limits  at  tg  as  shown  in  Figure  III-3  and  to 
utilize  the  integral  equation  for  momentum  derived  in  Appendix  B.  Since 
there  are  no  external  forces  or  any  imposed  pressure  gradient,  the 
applicable  momentum  equations  reduce  to  the  following: 


where 


T = the  stress  tensor. 

n = the  unit  vector  perpendicular  to  subregion  surface 
as  shown  in  Figure  III-4. 

A = i 3/3x  + j 3/3y. 

p = density. 


For  this  problem,  the  stress  tensor  reduces  to  the  following 


Consequently,  the  term 


becomes 


where  for  simplicitity,  the  value  of  the  viscosity  p is  set  equal  to  1 


For  this  two  dimensional  problem,  the  second  term  on  the  right  side 
of  Equation  III-3  are  line  integrals.  Since  the  stress  tensor  reduces 
to  Equation  I I 1-4,  the  terms  in  question  are,  in  this  case,  independent 
of  the  y coordinate.  Therefore,  only  integrations  along  the  top  and 


bottom  of  the  subregions  are  needed. 

In  order  to  carry  out  the  calculations,  the  adjacent  cells  must  be 
inter-related  through  an  averaging  of  the  stress  tensor.  In  Figure 
I1I-5,  three  adjacent  cells  are  presented,  where  the  slopes  of  the 
velocities  are  indicated  by  and  (aj^)j_j.  The  procedure 

for  obtaining  the  average  consist  of  assuming  that  these  velocity  slopes 
exist  at  the  center  of  each  cell  and  then  fitting  them  with  a least 
square  cubic  equation  of  the  following  form: 

**  2 

a^  = y + C2  y + Cj  . (III-5) 


There  are  three  boundary  conditions  which  are  accounted  for  as  dia- 
gramed in  Figures  III-6A  to  III-6B.  After  the  momentum  distribution  at 
tg  + At  is  computed  for  all  of  the  non-boundary  cells,  the  boundary 
conditions  are  imposed.  The  upper  plate  remains  motionless,  thus  the 
fluid  velocity  adjacent  to  the  plate  remains  zero.  Consequently,  for 
the  upper  boundary  cell,  the  momentum  distribution  is  set  equal  to  a 
linear  fit  between  zero  at  the  top  of  cell  CL-j.  and  the  momentum  value  at 
the  top  of  the  next  lower  cell  CLj_,.  The  lower  plate  is  forced  to 
retain  a constant  velocity  of  0.8,  thus  for  the  lower  boundary  cell,  the 
momentum  distribution  is  set  equal  to  a linear  fit  from  the  momentum 
value  at  the  bottom  of  the  next  higher  cell  constant  plate 

velocity  at  the  bottom  of  the  lower  boundary  cell  CLp.  In  addition,  as 
the  momentum  is  propagated  in  the  y direction,  there  will  be  some  cells 
with  no  momentum.  Therefore  a moving  boundary  condition  exists  until 
all  of  the  cells  between  the  plates  have  gained  some  momentum.  In  order 
to  handle  this  moving  boundary,  the  cell  with  no  momentum,  say  CLj  which 
is  adjacent  to  a cell  with  momentum,  say  CL._  .,  is  identiHed.  Then  a 
momentum  distribution  is  assigned  to  CLr  whicn  is  a linear  fit  between 
zero  momentum  at  the  top  of  cell  CL™  and  the  momentum  value  at  the  top 
of  cell  CLj  j. 

Figure  I I 1-7  presents  the  momentvim  distributions  as  a function  of  y 
for  various  instants  of  time.  As  the  graph  shows,  momentum  was  pro- 
gressively transported  in  the  y direction  until  the  fluid  in  all  of  the 
cells  obtained  values  of  momentum  above  zero.  Then  the  momentum  in- 
creased in  all  of  the  cells  until  the  steady  distribution,  represented 
by  the  dashed  line,  was  obtained.  The  x's  on  the  graph  represent  the 
results  of  the  last  cycle  of  computation.  The  data  presented  in  Table 
III-l  represents  the  momentum  computed  at  the  top  and  bottom  of  the  cells 
where  the  values  were  determined  from  the  momentum  distribution.  The 
data  indicates  the  closeness  of  the  values  of  momentum  between  adjacent 
cells  and  suggests  that  a high  degree  of  continuity  between  cells  exist 
using  the  proposed  method. 
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III-6C  Lower  Plate  Boundary  Condition. 
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Figure  I I 1-7 


Momentum 


Momentum  Distributions  As  A Function  Of  Time  For  The 
Sliding  Parallel  Flat  Semi-infinite  Plate  Problem. 
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Table  III-l.  Momentum  Values  At  The  Top  and  Bottom  of  Cells  For 
The  Sliding  Plate  Problem.* 


Momentum  Values 


Cell 

Bottom 

Top 

1 

0.8 

0.729246870 

2 

0.729246870 

0.659359231 

3 

0.659359241 

0.591170446 

4 

0.591170426 

0.425451446 

5 

0.525451478 

0.462884307 

6 

0.462884251 

0.404040364 

7 

0.404040364 

0.349365223 

8 

0.349365239 

0.299169344 

9 

0.299169299 

0.253625760 

10 

0.253625788 

0.212773912 

11 

0.212773847 

0.176527839 

12 

0.176527922 

0.144689503 

13 

0.144689462 

0.116963056 

14 

0.116963063 

0.092972181 

15 

0.092972282 

0.072275841 

16 

0.072275821 

0.054383357 

17 

0.054383268 

0.038768292 

18 

0.038768374 

0.024879464 

19 

0.024879485 

0.012149046 

20 

0.012149046 

0.0 

★ 

Ax  = 1;  Ay  = 1;  Atime  = 0.1;  Time  = 40 
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IV.  A CONCEPTUAL  DESCRIPTION  OF  THE  CODE  CONSTRUCTION 


I'he  most  efficient  device  for  organizing  a conceptual  description 
of  a code  structure  is  the  flowchart  as  presented  in  Figures  IV-1  to 
lV-(5.  This  flowchart  is  a preliminary  effort  and  is  not  intended  to 
accomplish  more  than  a demonstration  of  an  overall  approach. 

That  part  of  the  flowchart  in  Figure  lV-1  symbolizes  the  modeling 
of  a real  problem  in  terms  of  the  actual  data  required  to  compute  a 
solution.  In  general,  such  a model  would  be  based  on  the  various  physi- 
cal characteristics  of  the  real  target  and  the  offensive  mechanism  with 
which  the  target  will  be  engaged.  In  particular,  the  geometry  of  the 
target,  the  properties  of  the  materials  of  which  the  target  is  constructed, 
and  any  specific  initial  conditions  must  be  defined.  In  addition,  the 
offensive  mechanism  must  be  characterized  and  described  in  an  idealized 
manner  appropriate  for  computer  computations.  Based  on  this  information, 
the  actual  input  data  for  a calculation  is  to  be  formalized. 

The  input  model  would  consist  primarily  of  basic  assumptions, 
boundary  conditions,  and  other  basic  data  associated  with  the  target 
materials.  The  basic  assumptions  deal  with  items  such  as  the  number  of 
dimensions,  compressibility,  the  number  of  kinds  of  species,  external 
forces,  and  etc.  The  boundary  conditions  would  vary  from  problem  to 
problem  and  accounting  for  them  will  constitute  one  of  the  more  difficult 
tasks.  Some  of  the  other  basic  data  referred  to  above  could  be  vapor- 
ization rates,  reaction  rates,  an  ignition  criteria  and  etc,  which  would 
be  needed  if  the  problem  required  accounting  for  combustion. 

Once  the  input  model  has  been  read  into  the  computer,  the  code 
should  execute  those  subroutines  designed  to  compute  moments  of  the 
physical  parameters  for  all  of  the  cells  and  then  sum  them  in  the  appro- 
priate fixed  cells.  These  subroutines  are  indicated  in  that  part  of  the 
flow  chart  shown  in  Figures  IV-2  and  IV-3;  beginning  at  SI  and  ending 
with  S.3. 

Assuming  that  the  cells  are  numbered  consecutively  from  1 to  a 
maximum  number  In,ax>  code  would  naturally  begin  with  the  1=1  cell 
and  resolve  the  question  as  to  whether  it  is  a boundary  cell.  If  the 
cell  is  one  defining  a boundary,  the  code  would  activate  that  subroutine 
(not  shown  on  the  flowchart)  which  accounts  for  the  appropriate  boundary 
condition  corresponding  to  that  cell,  and  then  proceed  to  the  next  I 
cell. 


Whenever  the  cell  does  not  constitute  one  defining  a boundary  con- 
dition, the  code  would  execute  a subroutine  entitled  "Neighboring  Cells" 
in  the  flowchart.  This  subroutine  would  identify  those  cells  which  are 
adjacent  to  the  current  I cell  by  their  I number  and  these  numbers  would 
be  stored  for  future  use.  The  reason  this  information  will  be  needed  is 
that  the  Lagrangian  cell  will  deposit  its  "contents"  in  certain  ones  of 
these  neighboring  cells,  depending  on  the  direction  of  the  fluid  velocity. 


I’.AI.CIILA'I  IONS  FOR  MASS,  MOMENTUM,  AND  ENERGY  MOMENTS  ( CONTINUED  J 
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; Figure  lV-4  Flowchart  For  Code  Construction  (Continued). 
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CALCUI-ATIONS  FOR  NF.W  DISTRI BIH  IONS 
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Figure  IV-5  Flowchart  For  Code  Construction  (Continued) 
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In  the  case  of  two  dimensional  flow,  there  will  be  nine  neighboring 
cells  including  the  current  I cell. 

The  next  step  is  to  determine  if  there  are  more  than  one  kind  of 
species.  In  the  case  of  one  kind  of  species,  the  code  would  proceed  to 
S2  of  the  flowchart  as  shown  in  Figure  lV-3.  If  not,  then  the  code 
would  initiate  a reiteration  design  to  compute  the  moments  of  the  species 
concentration  for  each  i kind  of  species.  Each  kind  of  species  will 
have  associated  with  it  a velocity  distribution  vi.  Consequently,  a 
Lagrangian  cell  of  that  species  would  be  calculated  to  move  according 
to  v^  and  limits  of  subregions  SR^  at  t^  would  be  determined  by  a sub- 
routine entitled  "Limits"  in  the  flowchart.  In  addition,  subroutine 

"Limits"  would  identify  the  neighboring  cells  in  which  the  moments  of  the 
subregions  are  to  be  summed.  Then  taking  each  subregion  in  turn,  the 
moments  for  the  species  concentration  would  be  calculated  and  summed. 

In  thi^  procedure  some  of  the  subroutines  would  be  bypassed  if  Wj  = 0 
or  if  fj^  = 0.  After  all  of  the  kinds  of  species  have  been  considered, 
the  code  would  move  on  to  consider  the  remaining  physical  parameters, 
where  the  flowchart  is  marked  with  the  symbol  S2. 

In  considering  the  computation  of  moments  for  the  remaining  physical 
parameters,  in  Figure  IV-3,  it  is  noted  that  subroutine  "Limits"  would 
again  be  called  and,  based  on  the  mass  average  velocity,  the  limits  of 
subregions  at  tg  would  be  determined.  Following  the  flowchart,  it  can 
be  seen  that  if  the  fluid  is  assumed  to  be  incompressible,  the  integral 
for  computing  the  mass  moments  would  be  bypassed.  Then  the  terms  for 
calculating  the  moments  for  momentum  and  energy  would  be  evaluated. 
However  certain  terms  may  be  bypassed  if  the  viscosity  is  zero,  the 
heat  flux  is  zero,  the  number  of  kinds  of  species  is  greater  than  one, 
or  if  exterior  forces  are  absent.  Once  all  of  the  subregions  have  been 
accounted  for,  the  code  would  recycle  until  all  of  the  cells  have  been 
processed.  At  the  conclusion  of  these  calculations,  corresponding  to 
each  of  the  I cells,  there  will  be  sums  of  the  various  moments  stored. 

The  new  distributions  of  the  physical  quantities  at  tg  + At  would 
be  determined  by  the  operations  presented  in  Figure  IV-5.  As  before  in 
the  case  of  boundary  cells,  special  subroutines  would  be  called  to  im- 
pose the  appropriate  boundary  conditions.  If  the  size  of  the  cells 
varies  from  one  cell  to  another,  then  for  each  cycle  of  computation  the 
Geometric  Moments  would  have  to  be  computed  for  the  different  size  cells 
or  a table  of  these  values  required.  If  the  cell  sizes  vary  in  the 
computation,  then  appropriate  adjustments  in  the  table  would  be  needed 
to  reflect  these  changes.  At  the  completion  of  these  calculations,  the 
new  distributions  of  the  physical  parameters  at  tg  + At  would  be  avail- 
able for  the  next  cycle  of  computation. 

In  Figure  IV-6,  some  of  the  operations  which  would  be  needed  to 
complete  the  calculation  are  presented.  Firstly,  certain  physical 
quantities  may  be  required  which  depend  on  those  already  calculated. 

These  could  include  fluid  temperature  and  pressure.  If  more  than  one 
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kilui  of  species  is  possible,  then  based  on  the  ignition  criteria,  it 
would  bo  determined  if  ignition  had  occurred.  If  so,  then  a new  value  of  w. 

(the  mass  rate  of  production  term)  based  on  the  existing  temperature 
and  species  concentrations  would  be  evaluated  for  use  in  the  next 
cycle.  Afterwards  a new  species  velocity,  for  each  species,  would  be 
determined  from  which  a new  mass  average  velocity  could  be  derived.  At 
this  point,  the  code  would  be  prepared  to  print  out  data,  initiate  the 
next  cycle  or  stop. 

Ihe  flowchart  does  not  indicate  numerous  operations  and  subroutines 
which  the  final  version  of  the  code  will  need.  In  some  cases  it  will 
be  desirable  to  change  the  cell  sizes  by  doubling  or  halving  some  or  all 
of  them.  A difficulty  will  be  the  programing  of  the  code  such  that 
arbitrary  coordinates  systems  can  be  handled  (cartesian,  cylindrical, 
spherical,  etc).  Also  certain  problems  must  be  resolved,  such  as  ensur- 
ing that  distributions  of  physical  quantities  are  continuous  across  cell 
boundaries  and  that  non-negative  quantities  never  end  up  with  negative 
values.  In  order  to  place  these  and  possibly  other  problem  areas  in 
perspective  relative  to  the  entire  code,  the  flowchart  needs  to  be 
expanded.  However,  the  current  version  provides  sufficient  descriptive 
information  to  adequately  portray  the  basic  approach  for  organizing  the 
code  structure. 


V.  CONCLUSIONS 

The  objective  of  this  report  is  to  document  the  essential  and  basic 
elements  of  the  proposed  mathematical  and  physical  scheme  for  simulating 
fluid  flow.  The  simple  problem  solved  in  Section  III  provides  verifi- 
cation to  the  extent  that  it  demonstrates  that  the  approach  is  feasible. 
The  flowchart  and  its  associated  description  in  Section  IV  provides  a 
general  picture  as  to  a possible  procedure  for  implementing  the  proposed 
method  in  terms  of  a computer  code.  There  is  no  attempt  to  anticipate 
or  categorize  all  of  the  difficult  areas  or  to  explain  how  they  might  be 
overcome.  That  would  tend  to  confuse  and  undermine  the  primary  objective 
of  the  report. 

The  proposed  method  has  two  important  characteristics.  One  is  that 
it  involves  integral  equations,  which  means  that  singular  points  will 
not  be  a problem  and  no  instabilities  will  occur  for  that  reason.  The 
other  is  that  the  various  physical  parameters  will  be  represented  by 
variable  distributions  across  the  cells,  which  should  provide  increase 
accuracy  in  contrast  to  using  constant  distributions  for  the  same  cell 
size. 

At  this  stage,  it  is  impossible  to  make  a reasonable  comparison  of 
this  method  with  existing  hydrodynamic  codes  which  utilize  other  methods. 
This  can  be  accomplished  only  when  a code  has  been  written  and 
computational  results  compared.  There  is  no  way  to  estimate  the  computer 
time  required  to  compute  a typical  problem.  That  also  must  be  held  in 
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abeyance  until  problems  are  run.  To  this  end,  progress  in  code  construc- 
tion will  be  accomplished  as  manpower  and  funds  become  available. 


THE  DETERMINATION  OF  DISTRIBUTIONS  OF  FUNCTIONS  BY 
THE  MOMENTS  METHOD 


\ 

• s 


A part  of  the  mathematical  basis  of  the  theory  is  the  so  called 
Moments  Method  for  the  determination  of  an  average  distribution  of  a 
variable  over  a region  if  certain  moments  of  the  variable  over  sub- 
regions  are  knovm.  The  method  may  be  derived  directly  from  the  concept 
of  a weighted  arithmetic  mean.  In  general,  if  it  is  assumed  that  there 

are  available  a^,  02,  Oj  non-negative  numbers  (weights),  not  all 

zero,  the  weighted  arithmetic  mean  f of  fj,  f2,  fj  is  defined  by 

the  following  formula: 


f 


(A-1) 


IVhen  the  weights  are  all  equal.  Equation  A-1  reduces  to  the  expression 
for  computing  the  ordinary  arithmetic  means. 

This  may  be  extended  to  the  weighting  of  a function  over  an  interval 
where  the  weighting  function  varies  continuously.  This  extention  yields 
the  following  equation: 


/ a(x)  f(x)  dx 

f(x)  = ^ 

/ a(x)  dx 


(A-2) 


where  a(x)  is  a non-negative  weighting  function  whose  integral  is  not 
zero. 

In  Figure  A-1,  a function  a(x),  whose  first  derivative  is  discon- 
tinuous at  certain  values  of  x is  plotted.  Between  each  pair  of  first 
derivative  discontinuity  values,  the  function  f(x)  and  the  weight  func- 
tion a(x)  are  assumed  to  be  known.  Utilizing  this  information  and 
Equation  A-2,  the  following  expression  can  be  obtained: 


Also  indicated  in  Figure  A-1,  by  the  dashed  line,  is  the  average 
weight  function  for  the  entire  range  of  x.  The  integral  of  the  average 
weight  function  over  the  entire  range  of  x is  equivalent  to  the  sum  of 
the  integrals  of  the  individual  weight  functions  over  their  correspond- 
ing intervals  of  x.  That  is 


we  may  write 


where  a ^ and  f are  known.  If  f is  set  equal  to  x = 1,  then  a is  a 
constant  as  indicated  in  Figure  A-1.  However,  f may  be  any  power  of  x 
and  for  that  reason  o can  be  evaluated  as  a variable  function.  In  two 
or  three  dimensions,  f may  be  a function  of  the  other  spatial  parameters 
as  well. 


For  the  purpose  of  demonstration,  we  may  utilize  Figure  A-2,  where 
a two  dimensional  spatial  region  is  defined  and  which  is  partitioned 
into  five  subregions  SRj.  Also,  we  may  assume  that  corresponding  to 
each  subregion  a known  distribution  of  a physical  parameter  a(x,y')  exist 


and  dA  = dx  dy. 

The  integrals  on  the  left  of  Equations  A-9  can  be  evaluated  once 
the  limits  of  the  subregions  have  been  determined.  By  summing  the  values 
of  those  integrals  and  substituting  from  Equations  A-8,  the  following 
form  of  the  equations  is  obtained: 
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r 


k 


i,  1 


J SR 


/ “l  “ ■ "l  /l-  “ • >>3  /dA 

SRj  TC  TC  TC 


(A-10) 


^ J X a , dA  = bj  J'x^  dA  + b2  y dA  + bj  J’x  dA 

J SRj  TC  TC  TC 

E />"  “j  " “ '•i  y dA  + b^  dA  . b3  dA 

J SR,  TC  TC  TC 


The  integrals  on  the  right  side  of  Equations  A-10  depend  entirely 
on  tlie  spatial  parameters,  so  for  convenience,  these  integrals  are  re- 
ferred to  as  Geometric  Moments  and  are  integrated  over  the  entire  spatial 
region.  Since  there  are  three  equations  and  three  unknowns,  the 
coefficients  bj,  b2,  and  b,  of  a can  be  evaluated  by  solving  the  equa- 
tions simultaneously.  If  higher  order  expressions  for  a are  desired, 
the  additional  coefficients  required  can  be  evaluated  by  increasing  the 
value  of  the  k's  in  Equation  A-7  and  thereby  increasing  the  number  of 
Equations  in  A-10  accordingly. 
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APPENDIX  B 

THE  DERIVATION  OF  THE  INTEGRAL  EQUATIONS  FOR  CALCULATING 
MOMENTS  IN  FLUID  FLOW 


In  deriving  the  integral  equations  for  computing  moments  of  physical 
quantities  over  a spatial  domain,  where  the  results  include  changes 
over  a time  step  At,  it  is  sufficient  to  begin  with  the  basic  equations 
of  change,  written  as  follows: 


Species  Concentration: 


Mass : 


Momentum: 


3n.  w. 

+ V • n.v.  = — 

at  1 1 m^^ 


■^  + V • pv  = 0 

at 


+ 7 . 


7 • (pv)  v=-7P-  (7*t)+E  n.m.f. 


Total  Energy: 


3E  r,  -►cT 
~ + 7 . vE 


Kinetic  Energy: 


3E  r, 

. V . E V . 


- 7 • (t  • - (7  • Pv)  - (7  • q) 


+ E n.m.?.  • 

i 1 1 1 1 


4 -► 

7 • (t  • v)  + (t  : 7 v)  + P(7  • v) 


- 7 • Pv  + E n.m. f . • V. 

4 111  1 


Internal  Energy: 


— + 7 . E V = 


-7‘q-T:7v-P(7*v) 


(B-la) 


(B-lb) 


(B-lc) 


(B-ld) 


(B-le) 


(B-lf) 
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The  definitions  of  the  physical  quantities  used  in  Equations  B-1 
are  listed  as  follows: 

n^  - the  number  desity  of  the  species. 

- the  mass  rate  of  production  of  the  i^^  species 
through  chemical  reactions. 

th 

Vj  - the  average  velocity  of  the  i species. 

m^  - the  mass  of  the  i^*'  species. 

p - the  mass  per  unit  volume  (density) . 

V - the  mass  average  velocity  of  the  mixture. 

P - the  pressure. 

T - the  stress  tensor. 

th 

f.  - the  external  force  per  unit  mass  on  the  i 
species . 

T 

E - the  total  energy, 
q - the  energy  flux. 

E*^  - the  kinetic  energy  (1/2  p v^) . 

E^  - the  internal  energy. 


The  conservation  laws,  in  integral  form,  are  to  be  applied  to  fluid 
volumes  which  move  through  a fixed  coordinate  system.  In  addition, 
weight  functions  are  utilized  and  these  are  assumed  to  be  products  of 
powers  of  x,  y,  and  z at  t = tg  + At  and  in  addition  satisfy  the  follow- 
ing differential  equation. 


V . V f = 0 


(B-2) 


In  deriving  the  necessary  form  of  the  equations,  the  function  f is 
initially  multiplied  through  Equations  B-1.  Substitution  and  rearrange- 
ment of  the  terms  yields  the  following  form: 

Species  Concentration: 


3fn.  - w. 

V + f (V  • (n.  V.)  = n.  ^ + f — 
9t  ' ^11-'  1 at  m. 

1 


(B-3a) 


Moment  urn 


Total  Energy 


Kinetic  Energy 


Internal  Energy 


The  following  identities  can  be  used  for  transforming  the  above 
equations  to  the  desired  form: 


and 


(B-4b) 


V • UVv  = V(V  U • v)  ♦ U(7  • v)  ♦ UV(V  • v) 


U V V • aT  = (V  • UVv)  - (V  V U • V)  - (UV  V • V) 
V • UVv  = UV(V  • V)  + (^  • V UV) 


where  U and  V are  arbitrary  functions. 


(B-4c) 


Substitution  of  Equations  B-4  into  Equations  B-3  and  after  rearrange 
ment  of  terms,  the  following  equations  are  obtained: 

Species  Concentration: 

3fm.  -o  w 

+ V.  • V fn.  + fn.  (V  . V.)  = n.  + V.  • Vf)  + f — (B-5a) 


Mass : 


1^  + V • Vfp  + fp  (V  • v)  = pC^I^  + V • Vf) 


(B-Sb) 


Momentum : 


+ V • Vfpv  + f pv  CV  • v)  = p V ^ 


- f (V  • t)  - f VP  + f E n^  m^  r 


(B-5c) 


Total  Energy: 

^ + V • Vf  E^  + f E^  (V  • v)  = E^  (|^  + V • Vf)  - f (V  • (?  • v)) 


(B-5d] 


- f(V  • q)  - f(V  • P aT)  + f E n,  m.  f.  • v. 


Kinetic  Energy: 


+ V • (V  . v)  - (||  * V • Vf) 


f V • (?  • V)  + f (f  : V ^)  + f P (V  • v) 
f (V  • P^)  + f E ra^  fj  • 


Internal  Energy; 


— + V • VfE^  (V  • v)  = tlf  "■  V • 
3 t 


- f (V  • q)  - f (t  : 7^  - f P (V  • 


(B-5e] 


(B-5f) 


The  substantial  derivative  (also  referred  to  as  the  derivative 
following  the  motion)  D/DT  is  defined  as  the  change  with  respect  to 
time  plus  the  change  with  respect  to  position.  That  is: 


Huvi  ^ i(yvi  ^ ^ uv 

Dt  3t 


(B-b) 


Substitution  of  the  substantial  derivative  into  Equations  B-5  and 
then  integrating  the  terms  over  a spatial  domain  R(t),  which  has  a 
bounding  surface  R*(t),  yields  the  following  set  of  integral  equations: 


Species  Concentration: 

f (V  . 7p1dR.  /"iflf  * 

Kt)*-  J Rtt)  L J (B-7a) 

+ /*  f — dR 

J 
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Mass: 


/[ 

R(t) 


Pfp 

bt 


+ f p (V  • v) 


1 <"*'  / » [if  * • ">1 

R(t)  •- 


+ (V  • Vf)  dR  (B-7b) 


Momentum: 

Dfpv 
Dt 


R(t)  '■ 


+ f p V (V  • v)ldR  = 


) dR  = y p vlll  + • Vf)ldR 

R(t)  L J 


y f (\/  • ?)  dR  - y f V P dR 
R(t)  R(t)  R(t) 


y f V P dR  + y f E m^ 


(3-7c) 


dR 


Total  Energy: 


yfe. fE-t (, . ?)idR.  y .7 • vfidR 

R(t)  ^ R(t)  L 


ft7  ■ 


(t  • V)  dR  - y f (V  • q)  dR 
R(t)  R(t) 


/ 


(B-7d) 


/ 


f (V  • P v)  dR  + / f E n.  m.  f.  • V. 

-/■111! 

Rft)  R(t) 


dR 


Kinetic  Energy: 


/[ 

R(t) 


^ - eeSv 


/ 


• V)  j dR  = y + (V  * Vfj 


* / 


dR 


f [V  • (t  • v)]  dR  + J f (t  : V v)  dR 
R(t)  R(t) 


(B-7e) 


(continued) 
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*/ 


f P (V  • v)  dR  - J f (V  • P v)  dR 
R(t)  R(t) 


/ 


(B-  7c 
cont) 


/f  E n.  m.  • V. 
Jill  1 

R(t) 


dR 


Internal  Energy: 


/fe-  * f e'  (V.  7ldR  - / e'  [if  ♦ ? • 7fj  dR 

RCt)*-  -I  R(t)  *-  J 


/ 


f (V  • q)  dR  - y f (T  : V v)  dR 
R(t)  R(t) 


/f  (?  : 


/ 

R(t) 


f P (V  • v)  dR 


In  order  to  transform  Equations  B-7,  such  that  appropriate  integral 
equations  can  be  obtained,  it  is  necessary  to  recall  the  Convection 
Theorm^,  given  as  follows: 


^ y*  U dR  = yfe  + U (V  • v)!  dR  (B-8) 

R(t)  RCt)*- 

^ dR  = ^ JudR*  j U (V  • it)  dR*  (B-&1 

R(t)  R(t)  R(t) 


^Meyep,  Richard  E.,  "Introduotion  to  Mathematioal  Fluid  Dynamias," 
Wile  Interaaienoe,  New  York^  N.  Y. , 1971. 
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where  n is  a unit  vector  normal  to  the  surface  dR*.  For  our  purpose, 
the  second  term  of  Equation  B-9  is  set  to  zero  since  no  fluid  is  to  be 
permitted  to  cross  the  boundary  of  the  cell.  The  application  of  this 
condition  plus  the  imposition  of  the  restriction  as  defined  by  Equation 
B-2  to  Equations  B-7  yields  the  following  set  of  equations: 


Species  Concentration; 


Mass: 


^ f 

dt  J 


f p dR 


Momentum: 


(B-lOa) 


(B-lOb) 


^ ^ f pv  dR  = - ^ f (V  • ?)  dR  - ^ f 7P  dR 

R(t)  R(t)  R(t) 


r 

+ I f.  Z n.  m.  f. 
J 1 j 1 1 1 


(B-lOc) 


Total  Energy: 


i ft  f tv.  (f. 


v)  dR  - / f fv  • q)  dR 

R(t) 


y f (7  • Pv)  dR  + J 
R(t)  R(t) 


/f  Z n.  m.  • V. 

1 1 1 1 


(B-lOd) 


Kinetic  Energy: 


J"  f dR  = - y f 7 . (?  . v)  dR  + /*  f (? : 

R(t)  R(t) 

+ y f P (7  • V)  dR  - y f (7  • Pv)  dR 
R(t)  R(t) 

- - 

+ / f £ n.  m.  f.  • V.  dR 

J Jill  1 


Internal  Energy: 


7v)  dR 


(B-lOe) 


^ f 

dt  J 


f E^  dR  = - / f (7  • q)  dR  - / f (f  : 7^)  dR 


f P (7  • dR 


(B-lOf) 


Approximate  solutions  can  be  obtained  by  replacing  Equations  B-10 
with  the  following  equations: 

Species  Concentration: 

y f n.  dR  - / f n.  dR  + At  /"f  — dR  (B-lla) 

R(t  + At)  R(t J R(t  ) 

^ 5 S 


J f p dR  = J' f p dR 

R(t  ♦ At)  R(t J 

^ s 


(B-llbJ 


Moment  lun 


i'otal  Energy 


Kinetic  Energy 


Internal  Energy: 

y*  f E^  dR  = y*f 


/ 


dR  - At  / f (^  • q)  dR 


R(tg  ^ At) 


RCt,) 


Ht^) 


J 


-At  / f (t  : V v)  dR  - At  / f P (V  • v)  dR 


RU^) 


/■ 


(B-llf) 


The  following  expressions  are  required  to  obtain  the  desired  form 
of  the  equations: 


/ 


V • U V dR  = y*  U (7  • v)  dR  + y*  7 U • v dR  (B-]2a) 
3)  RCt^)  R(tp 


and  the  Divergence  Theorem: 


(7  • v)  dR  = 


R(t3) 


R*(t3) 


y*  U V . n 


dR 


(B-12b) 


Application  of  Equation  B-12  to  Equations  B-11  yields  the  final  form 
of  the  integral  equations  written  as  follows: 


Species  Concentration: 

/ f "i  f "i 


dR  ♦ At 


RCtg  + At) 


R(t  J 


/w . 

f — 
m. 

RCt^) 


dR 


(B-13a) 


Mass : 


y*  f p dR » y* 


f p dR 


(B-13b) 


R(tg  + At) 


RCt^) 


Momentum : 

y*  fpvdR=  f fpvdR-At 


R(t^  + At) 


/ 


dR 


R*(tp 


At  y*  7f  • ? dR  - At  y*f  7P  dR  + At  J" f I 


(B-13c) 

dR 


R(t3) 


R(t^) 


RCt^) 
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^ s 


t *’• 


1 


Total  Bnergy: 

J f E*  dR  = / f E‘  dR  - At 


/ f (t  • v)  *11 


dR 


R(t^  + At) 


R(t,) 


R*(tg) 


+ At 


J" Vf  • (t  • 


/ 


v)  dR  - At  /f  q • n dR  + At  # Vf  • q dR 


R(tj.) 


R*(t^) 


R(tg)  (B-13d) 


- At  ^ f (V  • P v)  dR  + At  J" f Z 


dR 


RCt^) 


R(t3) 


kinetic  Energy: 

E*^  dR  = E*^  dR  - At  f i (t  • v)  dR* 


/ ^ 
r' + Ati 


R(t^  + At) 


R(t,) 


R*(t^) 


/ 


+ At  ! vf  • (t  • v)  dR  + At  / f (t  : V v)  dR 


RCt,) 


R(tp 


* it  f { i 


? (V  • v)  dR  - At 


RCt^) 


ft  ■ 


v)  dR 


R(t^) 


f - - 

+ At  / f Z n.  m.  f.  • V. 

y ill!  1 


dR 


RCt3) 


Internal  Energy: 


y f e'  dR  = /: 


f E^  dR  - At  / f q • n dR* 


R(tg  + At)  R(tg) 


R*(t^) 


+ At  ^ Vf  • q dR  - At 


q dR  - At  I f (t  : V v)  dR 


RCt^) 


R(t3) 


/ 


-At  / f P {V  • v)  dR  . 


R(t3) 


(B-13e) 


(B-13f) 
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These  integral  equations  are  sufficient  to  utilize  for  the  computa- 
tion of  approximate  values  of  moments  of  physical  quantities  related  to 
fluid  flow  over  a spatial  domain  in  a time  increment  At. 
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APPENCIX  C 

DEVELOPMENT  OF  EXPRESSIONS  FOR  CALCULATING  DENSITY, 

MOMENTUM  AND  KINETIC  ENERGY  DISTRIBUTION  AT 
tg  + At  BASED  ON  THE  MOMENTS 

The  form  of  the  expression  for  calculating  a distribution  of  a 
physical  parameter  over  a spatial  region  is  arbitrary  provided  the  values 
of  the  required  moments  are  available.  Since  the  limits  of  the  sub- 
regions  are  determined  by  the  velocity  distributions  and  it  is  desirable 
to  simplify  the  integrations,  the  velocity  distribution  are  assumed  to 
be  linear  in  form,  as  follows: 

v^  = i (A(l,l)  y + A(2,l)  z + A(3.1)) 


Vy  = j (A(l,2)  z + A(2,2)  X + A(3,2)) 


(C-1) 


where 


v^  = k (A(l,3)  X + A(2,3)  y + A(3,3)) 


^x’  ~ velocity  components 


i , j , k = unit  orthogonal  vectors 

A(i,j)  = coefficients 

X,  y,  z = spatial  parameters  in  fixed 
coordinates . 


The  coefficients  of  Ej^pressions  C-1  cannot  be  evaluated  directly, 
but  must  be  evaluated  from  a knowledge  of  the  density  distribution  and 
the  momentum  or  the  kinetic  energy  distribution.  Consequently,  the 
density  is  required  to  be  a constant  across  the  cell  and  the  momentum 
distribution  is  a linear  function  similiar  to  Expressions  C-1.  The 
form  of  the  kinetic  energy  distribution  reflects  the  velocity  squared 
term.  To  assist  in  writing  down  the  equations  the  symbol  G,  which 
represents  the  Geometric  Moments  as  presented  in  Table  C-1,  are 
utilized. 


Fable  C- 

-1.  The 

Integrals  : 

for  Calculating 

Geometric 

Moments* 

G(l, 

,1) 

f 0 0 

= Jy  z 

dRi 

G(1 

.2) 

f 0 0 

= /z  X 

dR^ 

G(1 

.3) 

r 0 0 

= Jx  y 

dR3 

^/y  z 

dRi 

G(2 

.2) 

= /z  X 

dR2 

G(2 

,3) 

f 1 0 
= Jx  y 

‘*«3 

C.(3, 

,n 

= /yV 

dR^ 

G(3 

,2) 

r 2 0 

= JZ  X 

dR2 

G(3 

,3) 

r 2 0 

= Jx  y 

dR3 

G14, 

,1) 

= /yV 

dRi 

G(4 

.2) 

/•  3 0 
= jz  X 

dR2 

G(4 

.3) 

f 3 0 
= Jx  y 

dR3 

G(5, 

,1) 

= /yV 

dRi 

G(5 

.2) 

f 4 0 
= /z  X 

dR2 

G(5 

.3) 

f 4 0 
= Jx  y 

dR3 

G{6, 

,n 

dRj 

G(6 

.2) 

f 0 1 
= /z  X 

dR2 

G(6 

,3) 

r 0 1 
= Jx  y 

dR3 

G(7. 

,1) 

= /y^z^ 

dR^ 

G(7 

,2) 

r 1 1 

= /z  X 

dR2 

G(7 

,3) 

r 1 1 
= Jx  y 

dR3 

G(8. 

.1) 

= /y^z^ 

dRi 

G(8 

.2) 

r 2 1 

= /z  X 

dR2 

G(8 

.3) 

f 2 1 
= Jx  y 

dR3 

G(9, 

1) 

= jy^z^ 

dRj 

Gl9 

,2) 

r 3 1 

= Jz  X 

dR2 

G{9 

.3) 

f 3 1 
= Jx  y 

dR3 

G(10, 

1) 

dRj 

G(10 

,2) 

r 4 1 

- Jz  X 

dR2 

G(10 

.3) 

f 4 1 
= Jx  y 

dR3 

G(ll. 

,1) 

dRi 

G(ll 

,2) 

r 0 2 

= Jz  X 

dR2 

G(ll 

.3) 

r 0 2 

= Jx  y 

dR3 

G(12, 

.1) 

= ^^z^ 

dRi 

G(12 

.2) 

f 1 2 
= Jz  X 

dR2 

G(12 

.3) 

r 1 2 

= Jx  y 

dR3 

G(13, 

.n 

= jy^z^ 

dRi 

G(13 

.2) 

f 2 2 
= Jz  X 

dR2 

G(13 

.3) 

r 2 2 
= Jx  y 

dR3 

G(14. 

,1) 

=^v 

dRi 

G(14 

.2) 

r 3 2 

= Jz  X 

dR2 

G(14 

,3) 

f 3 2 
= Jx  y 

dR3 

G(15, 

,1} 

r 4 2 
= Jy  z 

dRi 

G(15 

.2) 

f 4 2 
= Jz  X 

dR2 

G(15 

,3) 

f 4 2 
= /x  y 

dR3 

G(16, 

,1) 

C 0 3 

= Jy  z 

dRi 

G(16 

.2) 

r 0 3 
= Jz  X 

dR2 

G(16 

,3) 

f 0 3 
= /x  y 

dR3 

G(17, 

,1) 

t 1 3 
= jy  z 

dRj 

G(17 

.2) 

f 1 3 
- Jz  X 

dR2 

G(17 

,3) 

r 1 3 
= /x  y 

dR3 

G(18, 

.1) 

f 2 3 
= /y  z 

dRj 

G(18 

.2) 

f 2 3 
-Jz  X 

dR2 

G(18 

.3) 

r 2 3 
= /x  y 

‘^«3 

G(19. 

,1) 

,3  3 
= _^  z 

dR^ 

G(19 

.2) 

f 3 3 
= Jz  X 

dR2 

G(19 

.3) 

r 3 3 
= /x  y 

G(20. 

.1) 

=^v 

dRj 

G(20 

.2) 

f4  3 
- jz  X 

dR2 

G(20 

.3) 

f 4 3 
= Jx  y 

G(21, 

.1) 

f 0 4 
= Jy  z 

dRj 

G(21 

.2) 

f 0 4 
= Jz  X 

dR2 

G(21 

,3) 

r 0 4 
= Jx  y 

•k 

dR^  = dx  dy  dz;  dR^  = dy  dz  dx;  dR^  = dz  dx  dy 
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Multiplying  through  by  f yields  the  following  expressions: 


f P 
f p V 

y 

-► 

f P V 

z 

where 


B(i,j) 


T (B(l,l)  f y + B(2,l)  f z ♦ B(3,n  f) 

T (B(l,2)  f z ♦ B(2,2)  f X ♦ B(3,2)  f) 

k CB(1,3)  f X + B(2,3)  f y + B(3,3)  f) 

the  coefficients 
velocity  components 

unit  orthogonal  vectors. 


(C-6) 


These  equations  dictate_^that  f be  assigned  (l,y,z)  for  tl^e  i 
Component,  (l,z,x)  for  the  j Component,  and  (l,x,y)  for  the  k Component. 
In  addition,  the  o^de^  of  integration  is  dx  dy  dz,  dy  dz  dx,  and  dz  dx 
dy  for  components  i,  j,  and  k respectively.  Consequently,  the  expres- 
sions for  obtaining  the  values  of  the  coefficients  of  the  three  com- 
ponents of  momentvun  are  as  follows: 

i Component: 

M(l,l)  = B(l,l)  G(2,l)  + BC2,n  G(6,1)  ♦ B(3,l)  G(l,l) 

M(2,l)  = B(l,l)  G(3,l)  + B(2,l)  G(7.1)  + B(3,l)  G(2,l)  (C-7a) 


M(3,l)  = B(l,l)  G(7,l)  + B(2,l)  G(ll,l)  + B(3,l)  G(6.1) 


where 


M(l,l  = H / (P  v^)j  dRj 
J SRj 

M(2.i)  =x;  / y (P  Vj 

J SRj 

M(3,l)  =23  / 2 (P  Vj  dRj 
J SRj 

and  the  G’s  are  substituted  from  Table  C-1.  Naturally,  the  coefficients 
B(l,l),  B(2,l),  and  B(3,l)  are  to  be  evaluated  by  solving  Equations 
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C-6a  simultaneously.  The  symbol  dRj  is  equal  to  dx  dy  dz  as  before, 
j Component: 

M(l,2)  - B(l,2)  G(2,2)  ♦ B(2.2)  G(6,2)  ♦ B(3,2)  G(l,2) 

M(2.2)  « B{1,2)  G(3,2)  + B(2,2)  G(7,2)  + B(3,2)  G(2,2)  (C-8a) 

M(3,2)  = B(l,2)  G(7,2)  ♦ B(2,2)  G(ll,2)  + B(3,2)  G(6,2) 

where 

M(l,2)  J (p  Vy)j  dR2 

J SRj 

M(2.2)  = J]  y z (p  Vy)j  dR2  (C-8b) 

JSRj 

MC3.2)  = 2]  / X (p  Vy)j  dR2 
J SRj 

and  the  symbol  dR2  is  equal  to  dy  dz  dx. 
k Component: 

M(l,3)  = B(l,3)  G(2,3)  + B(2,3)  G(6,3)  + B(3,3)  G(l,3) 

M{2,3)  = B(l,3)  G(3,3)  + 8(2,3)  G(7,3)  + B(3,3)  G(2,3)  (C-9a) 

M(3,3)  = 8(1,3)  G(7,3)  + 8(2,3)  G(ll,3)  + 8(3,3)  G(3,3) 

where 

M(l,3)  = 2:  /(P  v^)j  dR3 
J SRj 

M(2,3)  = 23  / (P  ‘^‘^3 
J SRj 

m(3.3)  = e / y (p  dRs 

J SRj 
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where  the  symbol  dR^  is  equal  to  dz  dx  dy. 


The  Kinetic  Energy 

The  kinetic  energy  is  not  a vector  quantity.  However,  since  it  is 
directly  related  to  the  velocity  which  is  a vector  quantity_j^  w^  will 
refer  to  components  of  kinetic  energy  corresponding  to  the  i , j , and  k 
components  of  velocity.  The  fact  that  kinetic  energy  is  equal  to 
1/2  p is  reflected  in  the  six  term  expressions  which  follow: 

T Component: 

= C(l,l)y^  + C(2,l)z^  + C(3,l)  + CC4,l)yz  + C(5,l)  + C(6,l)z 
j Component: 

55  P Vy^  » CCl,2)z^  + C(2,2)x^  + C(3,2)  + C(4,2)zx  + C(5,2)z  + C(6,2)x 

^ (C-10) 

k Component: 

Js  P v,^  = C(l,3)x^  + C(2,3)y^  + C(3,3)  C(4.3)xy  + C(5,3)x  + C(6,3)y 

z 

where  the  C's  are  the  coefficients  to  be  evaluated. 

Multiplying  through  Equations  C-10  by  the  function  f yields  the 
following  expressions: 

i Component: 

fih  P v^^)  - C(l,l)  f y^  ♦ C(2.1)  f z^  ♦ C(3.1)  f 


+ C(4,l)  f yz  + C(5,l)  f y + C(6,l)  f Z 


j Component: 


P = C(l,2)  f z^  + C(2,2)  f x^  + C(3,2)  f 
+ C(4,2)  f zx  + C(5,2)  f z + C(6,2)  f x 


(C-11) 


k Component: 


P v,^)  - C(l,3)  f x^  + C(2,3)  f y^  + C(3,3)  f 
z 

♦ C(4,3)  f xy  + C(5,3)  f x + C(6,3)  f y. 
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(C-12b 
cont . ) 


j Component 


where 


i 


N(1.2)  = fff  (h  p Vy2)j  dR2 
J SRj 

N(2.2)  =X:  {JJ  z (’SP  Vy^dR^ 

J SRj 

N(3.2)  = 2:  WJ  " ^ V^J  ‘"‘'2 

J SRj 

N(4.2)  = SJSJ  ZX  ih  P '''y^)j  ‘^•^2 
J SRj 

N(5,2)  =2]  ih  P Vy^)j  dR2 

J SRj 

N(6.2)  = 2:  iif  x^  (Js  p dR^ 

J SRj 

and  dR2  = dy  dz  dx. 


(C-13b) 


k Component : 

N(l,3)  = C(l,3)  G(3,3)  + C(2,3)  G(11.3)  + C(3,3)  G(1.3) 

+ C(4,3)  G(7,3)  + C(5,3)  G(2,3)  + C(6,3)  G(6,3) 

N(2,3)  = C(l,3)  G(4,3)  ♦ C(2.3)  G(12,3)  + C(3,3)  G(2,3) 

+ C(4,3)  G(8,3)  + C(5,3)  G(3,3)  + C(6,3)  G(7,3) 

N(3,3)  = C(l*3)  G(8,3)  + C(2,3)  G(16,3)  + C(3,3)  G(6.3)  (C-14a) 

+ C(4,3)  G(12,3)  + C(5,3)  G(7,3)  + C(6,3)  G(ll,3) 

N(4,3)  = C(l,3)  G(14,3)  + C(2.3)  G(17,3)  + C(3,3)  G(7,3) 

+ C(4,3)  G(13,3)  + C(5,3)  G(8,3)  + C(6,3)  G(12,3) 

N(5,3)  = C(l,3)  G(5,3)  + C(2,3)  G(13,3)  + C(3,3)  G(3,3) 

+ C(4,3)  G(9,3)  + C(5.3)  G(4,3)  + C(6,3)  G(8,3) 

(Continued) 
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where 


The  relationships  between  the  coefficients  of  the  momentum  expres 
sions  and  the  kinetic  energy  expressions  are  listed  in  Table  C-2. 


Table  C-2. 


Relationships  Between  The  Coefficients  of 

Distributions  at  t * Lt 
s 


i Component 


The  Density  Coefficient 
A' 

j Component 


k Component 


A(l,l) 

A(2,l) 

A(3,l) 


The  Velocity  Coefficients 

A(l,2) 

A(2,2) 

A(3,2) 


A(l,3) 

A(2,3) 

A(3,3) 


The  Momentum  Coefficients 


B(l.l)  = a'A(1,1) 
B(2,l)  = A'A(2.1) 
B(3,l)  = A'AC3,1) 


B(1.2)  = A'A(1,2) 
B(2,2)  = A'A(2.2) 
B(3,2)  = a'A(3,2) 


B(1.3)  = A'A(1,3) 
B(2,3)  = A'A(2,3) 
B(3,3)  = A’A(3,3) 


The  Kinetic  Energy  Coefficients 


C(i.i) 

= 1/2  A'Ad.D" 

C(2.1) 

= 1/2  A'A(2,1)^ 

C(3.1) 

= 1/2  A’A(3,1)^ 

C(4.1) 

= a'a(1,1)A(3,1) 

C(5.1) 

= a'A(1,1)A(3,1) 

C(6,l) 

= A'A(2.1)A(3,1) 

C(l,2)  = 1/2  A A(l,2)‘‘ 
C(2,2)  = 1/2  a'A(2,2)^ 
C(3,2)  = 1/2  a'A(3,2)^ 
C(4.2)  = a'A(1,2)A(2.2) 
C(5,2)  = a'A(1.2)A(3,2) 
C(6,2)  = A'A(2.2)A(3.2) 


C(l,3)  = 1/2  A'A(1,3) 
C(2,3)  = 1/2  a'A(2,3)^ 
C(3,3)  = 1/2  a'A(3,3)^ 
C(4,3)  - a'a(1,3)A(2,3) 
C(5,3)  = a'a(1,3)A(3,3) 
C(6,3)  = A'A(2,3)A(3,3) 
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LIST  OF  SYMBOLS 


a constant  density  distribution 


coefficients  for  velocity  distribution  at  t ♦ At 


the  average  slope  of  v across  three  adjacent  cells 
a least  square  cubic  fit. 


coefficients  for  velocity  distributions 


coefficients  for  momentum  distributions 


and  b_  - coefficients  for  distribution  of  arbitrary  parameter  o 


coefficients  for  Kinetic  Energy  distribution  at  t + At 


and  c_  - coefficients  of  least  square  cubic  fit  of  a^  of 
adjacent  cells. 


fixed  coordinate  derivatives 


dR-,  dR,  - equals  dx  dy  dz,  dy  dz  dx,  and  dz  dx  dy  respectively 


derivatives  over  volumes 


derivatives  over  surfaces 


substantial  derivative 


total  energy 


internal  energy 


kinetic  energy 


function  to  be  weighted 


Geometric  Moments 


Geometric  Moments 


consecutive  numbering  of  computational  cells 


maximum  number  of  computational  cells 
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LIST  OP  SYMBOLS 


unit  directional  vector 


unit  directional  vector 


unit  directional  vector 


exponents  for  function  f 


moments  for  momentum 


moments  for  kinetic  energy 


energy  flux 


denotes  columns  of  cells  in  fixed  coordinates 


S7  - identifies  parts  of  flowchart 


limits  of  subregions  at  t 


SR(t  + At)  - limits  of  subregions  at  t ♦ At 


denotes  rows  of  cells  in  fixed  coordinates 


time  step. 

an  instant  of  time  corresponding  to  start  of  cycle  of 
computation. 

an  instant  of  time  corresponding  to  start  of  problem, 
implies  integration  over  total  cell, 
an  arbitrary  function, 
an  arbitrary  function. 


mass  average  velocity 


species  velocity 


tt  H 


LIST  OF  SYMBOLS 


w. 

1 


X.  y,  z 
7 

a/3t 

P 

ot 

A 

P 


- velocity  components. 

- mass  rate  of  production  of  the  i^^  species  through 

chemical  reactions. 

- coordinates  in  fixed  coordinate  system. 

^ ^ ^ 

- the  del  operator;  i 3/3x  + j 3/3y  + k 3/3z. 

- partial  derivative  with  respect  to  time. 

- fluid  viscosity. 

- an  arbitrary  physical  parameter. 

- delta. 

- density  distribution. 


- stress  tensor. 


SUBSCRIPTS 

i - species 

I - cells  i.i  computational  grid 

J - denotes  SR's  which  are  summed  in  a fixed  I cell  at 

t + At 
s 

K - denotes  SR's  of  Lagrangian  cell  at  t^ 

K - maximum  number  of  subregions  of  Lagrangian  cell  at  t 

indjc  s 

S - denotes  subregions  in  Section  III 

s - denotes  time  at  start  of  cycle  of  computation 

T - denotes  subregions  in  Section  III 

0 - denotes  instant  of  time  corresponding  to  the  start 

of  problem 
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